#ifndef _SVD_
#define _SVD_

#include "moremath.h"
#include <vl/VLd.h>

class SVD{
private:
	int svdcmp(Matd &u, Vecd &w, Matd &v);
	void svdksb(Matd& u, Vecd& w, Matd& v, Vecd &b, Vecd &x);
	void nullspace(Vecd& w, Matd& v, Matd &vNull);
	void cleandiag(Vecd &w);

public:
	Matd A;	// A = U*W*VT
	Matd U;
	Vecd W;
	Matd V;
	SVD(Matd &a){
		A=a;
		U=A;
		svdcmp(U,W,V);
		cleandiag(W);
	}
	~SVD(){}

	// simply return the svd
	inline void getSVD(Matd &u, Vecd &w, Matd &v){
		u=U; w=W; v=V;
	}
	// solve Ax=b
	inline void solve(Vecd &b, Vecd &x){
		svdksb(U, W, V, b, x);
	}
	// compue the null space matrix
	inline void getNullSpace(Matd &vNull){
		nullspace(W, V, vNull);
	};
};

#endif
